print('filling material components arrays')
from define_problem_space_parameters import material_types
import numpy as np
from calculate_domain_size import nx,ny,nz
from initialize_fdtd_material_grid import eps_r_x,material_3d_space,deps_x,tau_k_x,eps_r_y,deps_y,tau_k_y
from initialize_fdtd_material_grid import eps_r_z,deps_z,tau_k_z,sigma_e_x,sigma_e_y,sigma_e_z,mu_r_x,mu_r_y,mu_r_z
from initialize_fdtd_material_grid import sigma_m_x,sigma_m_y,sigma_m_z
from common import fillZero
# creating temporary 1D arrays for storing 
# parameter values of material types
print('filling material components arrays')
nmaterial_types = len(material_types)
t_eps_r =np.zeros((nmaterial_types))
t_mu_r =np.zeros((nmaterial_types))
t_sigma_e =np.zeros((nmaterial_types))
t_sigma_m =np.zeros((nmaterial_types))
t_deps =np.zeros((nmaterial_types))
t_tau_k =np.zeros((nmaterial_types))
#创建临时一维数组来存储媒质类型的参数值
for ind in np.arange(nmaterial_types):
    t_eps_r[ind]=material_types[ind]["eps_r"]
    t_mu_r[ind]=material_types[ind]["mu_r"]
    t_sigma_e[ind]=material_types[ind]["sigma_e"]
    t_sigma_m[ind]=material_types[ind]["sigma_m"]
    t_deps[ind]   = material_types[ind]["deps"]
    t_tau_k[ind]   = material_types[ind]["tau_k"]

#当t_mu_r和t_sigma_r等于0时，给他们赋一个很小的值来防止出现除以0的错误
def lowwer(u):
    u[ np.abs(u)<1e-20 ] = 1e-20
lowwer(t_mu_r)
lowwer(t_sigma_m)
lowwer(t_deps)
lowwer(t_tau_k)

print('Calculating eps_r_x')
#eps_r_x(i,j,k) is average of four cells  (i,j,k),(i,j-1,k), (i,j,k-1), (i,j-1,k-1) 
eps_r_x[:nx,1:ny,1:nz] = 0.25*(
    t_eps_r[material_3d_space[:nx,1:ny,1:nz]-1]
    +t_eps_r[material_3d_space[:nx,:ny-1,1:nz]-1]
    +t_eps_r[material_3d_space[:nx,1:ny,:nz-1]-1]
    +t_eps_r[material_3d_space[:nx,:ny-1,:nz-1]-1])

deps_x[:nx,1:ny,1:nz]=0.25*(t_deps[material_3d_space[:nx,1:ny,1:nz]-1]
                             +t_deps[material_3d_space[:nx,:ny-1,1:nz]-1]
                             +t_deps[material_3d_space[:nx,1:ny,0:nz-1]-1]
                             +t_deps[material_3d_space[:nx,:ny-1,:nz-1]-1])

tau_k_x[0:nx, 1:ny, 1:nz] = 0.25 * (
    t_tau_k[material_3d_space[0:nx, 1:ny, 1:nz]-1] +
    t_tau_k[material_3d_space[0:nx, 0:ny-1, 1:nz]-1] +
    t_tau_k[material_3d_space[0:nx, 1:ny, 0:nz-1]-1] +
    t_tau_k[material_3d_space[0:nx, 0:ny-1, 0:nz-1]-1]
)
# eps_r_x(1:nx,2:ny,2:nz) = t_eps_r[material_3d_space(1:nx,2:ny,2:nz))
# deps_x(1:nx,2:ny,2:nz)=t_deps[material_3d_space(1:nx,2:ny,2:nz))
# tau_k_x(1:nx,2:ny,2:nz)=t_tau_k[material_3d_space(1:nx,2:ny,2:nz))                        
print('Calculating eps_r_y')     
#eps_r_y(i,j,k) is average of four cells  (i,j,k),(i-1,j,k), (i,j,k-1), (i-1,j,k-1) 
eps_r_y[1:nx, 0:ny, 1:nz] = 0.25 * (
    t_eps_r[material_3d_space[1:nx, 0:ny, 1:nz]-1] +
    t_eps_r[material_3d_space[0:nx-1, 0:ny, 1:nz]-1] +
    t_eps_r[material_3d_space[1:nx, 0:ny, 0:nz-1]-1] +
    t_eps_r[material_3d_space[0:nx-1, 0:ny, 0:nz-1]-1]
)

deps_y[1:nx, 0:ny, 1:nz] = 0.25 * (
    t_deps[material_3d_space[1:nx, 0:ny, 1:nz]-1] +
    t_deps[material_3d_space[0:nx-1, 0:ny, 1:nz]-1] +
    t_deps[material_3d_space[1:nx, 0:ny, 0:nz-1]-1] +
    t_deps[material_3d_space[0:nx-1, 0:ny, 0:nz-1]-1]
)

tau_k_y[1:nx, 0:ny, 1:nz] = 0.25 * (
    t_tau_k[material_3d_space[1:nx, 0:ny, 1:nz]-1] +
    t_tau_k[material_3d_space[0:nx-1, 0:ny, 1:nz]-1] +
    t_tau_k[material_3d_space[1:nx, 0:ny, 0:nz-1]-1] +
    t_tau_k[material_3d_space[0:nx-1, 0:ny, 0:nz-1]-1]
)

# eps_r_y(2:nx,1:ny,2:nz) = t_eps_r[material_3d_space(2:nx,1:ny,2:nz))
# deps_y(2:nx,1:ny,2:nz)=t_deps[material_3d_space(2:nx,1:ny,2:nz))
# tau_k_y(2:nx,1:ny,2:nz)=t_tau_k[material_3d_space(2:nx,1:ny,2:nz))  
print('Calculating eps_r_z')     
#eps_r_z(i,j,k) is average of four cells  (i,j,k),(i-1,j,k), (i,j-1,k), (i-1,j-1,k) 
eps_r_z[1:nx, 1:ny, 0:nz] = 0.25 * (
    t_eps_r[material_3d_space[1:nx, 1:ny, 0:nz]-1] +
    t_eps_r[material_3d_space[0:nx-1, 1:ny, 0:nz]-1] +
    t_eps_r[material_3d_space[1:nx, 0:ny-1, 0:nz]-1] +
    t_eps_r[material_3d_space[0:nx-1, 0:ny-1, 0:nz]-1]
)

deps_z[1:nx, 1:ny, 0:nz] = 0.25 * (
    t_deps[material_3d_space[1:nx, 1:ny, 0:nz]-1] +
    t_deps[material_3d_space[0:nx-1, 1:ny, 0:nz]-1] +
    t_deps[material_3d_space[1:nx, 0:ny-1, 0:nz]-1] +
    t_deps[material_3d_space[0:nx-1, 0:ny-1, 0:nz]-1]
)

tau_k_z[1:nx, 1:ny, 0:nz] = 0.25 * (
    t_tau_k[material_3d_space[1:nx, 1:ny, 0:nz]-1] +
    t_tau_k[material_3d_space[0:nx-1, 1:ny, 0:nz]-1] +
    t_tau_k[material_3d_space[1:nx, 0:ny-1, 0:nz]-1] +
    t_tau_k[material_3d_space[0:nx-1, 0:ny-1, 0:nz]-1]
)
# eps_r_z(2:nx,2:ny,1:nz) = t_eps_r[material_3d_space(2:nx,2:ny,1:nz))
# deps_z(2:nx,2:ny,1:nz) = t_deps[material_3d_space(2:nx,2:ny,1:nz))   
# tau_k_z(2:nx,2:ny,1:nz) = tau_k_z(material_3d_space(2:nx,2:ny,1:nz))  
print('Calculating sigma_e_x')
# sigma_e_x(i,j,k) is average of four cells (i,j,k),(i,j-1,k), (i,j,k-1), (i,j-1,k-1) 

sigma_e_x[0:nx, 1:ny, 1:nz] = 0.25 * (
    t_sigma_e[material_3d_space[0:nx, 1:ny, 1:nz]-1] +
    t_sigma_e[material_3d_space[0:nx, 0:ny-1, 1:nz]-1] +
    t_sigma_e[material_3d_space[0:nx, 1:ny, 0:nz-1]-1] +
    t_sigma_e[material_3d_space[0:nx, 0:ny-1, 0:nz-1]-1]
)

print('Calculating sigma_e_y')
# sigma_e_y(i,j,k) is average of four cells
# (i,j,k),(i-1,j,k), (i,j,k-1), (i-1,j,k-1) 
sigma_e_y[1:nx, 0:ny, 1:nz] = 0.25 * (
    t_sigma_e[material_3d_space[1:nx, 0:ny, 1:nz]-1] +
    t_sigma_e[material_3d_space[0:nx-1, 0:ny, 1:nz]-1] +
    t_sigma_e[material_3d_space[1:nx, 0:ny, 0:nz-1]-1] +
    t_sigma_e[material_3d_space[0:nx-1, 0:ny, 0:nz-1]-1]
)
                    
print('Calculating sigma_e_z')
# sigma_e_z(i,j,k) is average of four cells
# (i,j,k),(i-1,j,k), (i,j-1,k), (i-1,j-1,k) 
sigma_e_z[1:nx, 1:ny, 0:nz] = 0.25 * (
    t_sigma_e[material_3d_space[1:nx, 1:ny, 0:nz]-1] +
    t_sigma_e[material_3d_space[0:nx-1, 1:ny, 0:nz]-1] +
    t_sigma_e[material_3d_space[1:nx, 0:ny-1, 0:nz]-1] +
    t_sigma_e[material_3d_space[0:nx-1, 0:ny-1, 0:nz]-1]
)

                    
print('Calculating mu_r_x')
# mu_r_x(i,j,k) is average of two cells (i,j,k),(i-1,j,k)
mu_r_x[1:nx, 0:ny, 0:nz] = 2 * (
    t_mu_r[material_3d_space[1:nx, 0:ny, 0:nz]-1] *
    t_mu_r[material_3d_space[0:nx-1, 0:ny, 0:nz]-1]
) / (
    t_mu_r[material_3d_space[1:nx, 0:ny, 0:nz]-1] +
    t_mu_r[material_3d_space[0:nx-1, 0:ny, 0:nz]-1]
)
                    
print('Calculating mu_r_y')
# mu_r_y(i,j,k) is average of two cells (i,j,k),(i,j-1,k)
mu_r_y[0:nx, 1:ny, 0:nz] = 2 * (
    t_mu_r[material_3d_space[0:nx, 1:ny, 0:nz]-1] *
    t_mu_r[material_3d_space[0:nx, 0:ny-1, 0:nz]-1]
) / (
    t_mu_r[material_3d_space[0:nx, 1:ny, 0:nz]-1] +
    t_mu_r[material_3d_space[0:nx, 0:ny-1, 0:nz]-1]
)
                    
print('Calculating mu_r_z')
# mu_r_z(i,j,k) is average of two cells (i,j,k),(i,j,k-1)
mu_r_z[0:nx, 0:ny, 1:nz] = 2 * (
    t_mu_r[material_3d_space[0:nx, 0:ny, 1:nz]-1] *
    t_mu_r[material_3d_space[0:nx, 0:ny, 0:nz-1]-1]
) / (
    t_mu_r[material_3d_space[0:nx, 0:ny, 1:nz]-1] +
    t_mu_r[material_3d_space[0:nx, 0:ny, 0:nz-1]-1]
)
                    
print('Calculating sigma_m_x')
# sigma_m_x(i,j,k) is average of two cells (i,j,k),(i-1,j,k)
sigma_m_x[1:nx, 0:ny, 0:nz] = 2 * (
    t_sigma_m[material_3d_space[1:nx, 0:ny, 0:nz]-1] *
    t_sigma_m[material_3d_space[0:nx-1, 0:ny, 0:nz]-1]
) / (
    t_sigma_m[material_3d_space[1:nx, 0:ny, 0:nz]-1] +
    t_sigma_m[material_3d_space[0:nx-1, 0:ny, 0:nz]-1]
)
                    
print('Calculating sigma_m_y')
# sigma_m_y(i,j,k) is average of two cells (i,j,k),(i,j-1,k)
sigma_m_y[0:nx, 1:ny, 0:nz] = 2 * (
    t_sigma_m[material_3d_space[0:nx, 1:ny, 0:nz]-1] *
    t_sigma_m[material_3d_space[0:nx, 0:ny-1, 0:nz]-1]
) / (
    t_sigma_m[material_3d_space[0:nx, 1:ny, 0:nz]-1] +
    t_sigma_m[material_3d_space[0:nx, 0:ny-1, 0:nz]-1]
)
                    
print('Calculating sigma_m_z')
# sigma_m_z(i,j,k) is average of two cells (i,j,k),(i,j,k-1)
sigma_m_z[0:nx, 0:ny, 1:nz] = 2 * (
    t_sigma_m[material_3d_space[0:nx, 0:ny, 1:nz]-1] * 
    t_sigma_m[material_3d_space[0:nx, 0:ny, 0:nz-1]-1]
) / (
    t_sigma_m[material_3d_space[0:nx, 0:ny, 1:nz]-1] + 
    t_sigma_m[material_3d_space[0:nx, 0:ny, 0:nz-1]-1]
)